Non-afRne displacements in crystalline solids in the harmonic limit 
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A systematic coarse graining of microscopic atomic displacements generates a local elastic de- 
formation tensor D as well as a positive definite scalar x measuring non-affinity, i.e. the extent to 
which the displacements are not representable as affine deformations of a reference crystal. We 
perform an exact calculation of the statistics of x and D and their spatial correlations for solids at 
low temperatures, within a harmonic approximation and in one and two dimensions. We obtain the 
joint distribution Pix, D) and the two point spatial correlation functions for x ^-nd D. We show 
that non-afHne and affine deformations are coupled even in a harmonic solid, with a strength that 
depends on the size of the coarse graining volume Q, and dimensionality. As a corollary to our 
work, we identify the field, /i^, conjugate to x and show that this field may be tuned to produce a 
transition to a state where the ensemble average, (x), and the correlation length of x diverge. Our 
work should be useful as a template for understanding non-affine displacements in realistic systems 
with or without disorder and as a means for developing computational tools for studying the effects 
of non-affine displacements in melting, plastic fiow and the glass transition. 



I. INTRODUCTION 

Understanding the mechanical response of soft and dis- 
ordered solids il^ such as polymer gels fabric [3], 
foams [4], colloids [5], granular matter [6] and glasses [3 
|8] is challenging because one is often lead to questions 
that lie on the boundaries of classical elasticity theory [S]- 
lllj . For example, under external stress, particles i within 
a solid undergo displacements = — away from 
some chosen reference configuration to their displaced 
positions r^. In a conventional homogeneous solid, such 
displacements are affine, in the sense that they can be 
expressed as u.^ = DR^, where D = K^^cr is the deforma- 
tion tensor related to the external stress a via the tensor 
of elastic constants K. This is not true if the solid is 
disordered at a microscopic level. 

One of the principal sources of non-affinity is a space 
(and possibly even time) dependent elastic constant [H]. 
The local environment in a disordered solid varies in 
space, depending crucially on local connectivity or co- 
ordination such that the local displacement may not 
be simply related to the applied stress a. Such non-affine 
displacements are present even at zero temperature, are 
material dependent, and vanish only for homogeneous 
crystalline media without defects. 

In this paper we explore another, perhaps complemen- 
tary, source of non-affinity, namely that which arises 
due to thermal fluctuations and coarse graining. Elas- 
tic properties of materials emerge upon coarse graining 
microscopic particle displacements [T71E5] over a coarse- 
graining volume r2. A systematic finite size scaling anal- 
ysis of the Q. dependent elastic constants then yields the 
material properties in the thermodynamic limit [211 122] ■ 
Such a coarse graining procedure has been used to ob- 



tain elastic constants of soft colloidal crystals from video 
microscopy [201 HIl H3] as well as in model solids [T71 [32] ■ 
For distances smaller than the size of fi, particle dis- 
placements, in d dimensions, are necessarily non-affine 
since the local distortion D of f2 is obtained by project- 
ing the displacements of all N particles in SI into the 
d X d-dimensional space of affine distortions which, in 
general, is smaller than the full iVc?-dimensional config- 
uration space available. The generation of non-affinity 
X, defined as the sum of squares of all the particle dis- 
placements which do not belong to the projected space 
of affine distortions, is therefore a necessary consequence 
of the coarse-graining procedure. Here we take a detailed 
look at this process and present an exact calculation for 
the probability distributions and correlation functions for 
X and D for harmonic solids in d = 1 and 2 in the canon- 
ical ensemble. Our work allows us to identify the field 
conjugate to x, viz. h^, and we show that by tuning this 
hy. one may enhance non-affine fluctuations and cause 
a transition. At this transition all the moments of the 
probability distribution of x diverge, thereby disordering 
the solid isothermally. 

There are several reasons why we believe that our work 
may be useful. Firstly, the harmonic crystal is often the 
starting point for more realistic calculations of the elas- 
tic properties of solids and constitutes an ideal system 
to which simulation and experimental results [23] can 
be compared in order to quantify purely anharmonic ef- 
fects. Secondly, a coarse-grained theory for the mechani- 
cal properties of soft solids should contain both the elastic 
and non-elastic fields D and x'- our work may provide a 
hint on how such a theory may be constructed. Thirdly, 
we believe that it may be possible to extend our calcu- 
lations to systems with isolated defects or randomness. 
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thus extending the analysis of Ref. [H] to non-zero tem- 
peratures. Finally, our calculations may be used to devise 
new simulational strategies for understanding the influ- 
ence of non-affine fluctuations on the mechanical proper- 
ties of both crystals and glasses and, perhaps, shed more 
light on the nature of the glass transition itself. 

The paper is organized as follows. In section [Tl] we 
set up the calculation and define Xi ^ and the coarse- 
graining process. In section [lIl] we present our calculation 
for the single point probability distributions for x &nd D 
for d = 1 harmonic chains and the d = 2 triangular har- 
monic lattices. Approximate experimental realisations of 
these systems correspond to mercury chain salts |13j and 
the spectrin network in red blood corpuscles [2] respec- 
tively. In section IV we evaluate the spatial correlation 
functions for D and x- This is followed by a calculation 
of linear response and identification of the non-affine field 
in section |V| Finally we discuss our results and conclude 



by giving indications of future directions in section VI 



II. COARSE GRAINING AND THE 
NON-AFFINE PARAMETER 

Consider a neighborhood, fi, in a d dimensional lattice 
consisting of N particles i arranged around the central 
particle within a cut-off distance Rn ■ Mostly we set Rq 
equal to the nearest neighbour distance so that con- 
tains all nearest neighbours of particle 0, but in section [VT| 
we also consider larger Rq . The zero temperature lattice 
positions that we choose as our reference are Ri and Rq 
and the fluctuating atom positions will be denoted r j and 
To [ISJ. Define the particle displacements = — R^, 
and Ai = Uj — Uo = r,; — ro — (R^ — Ro) as the displace- 
ment of particle i relative to particle 0. We will often use 
the Fourier transform of the particle displacement, Uq, 
which is defined such that the real-space displacements 
are Aj = Uj — Uq = Ivj^^ J (iqUq(e*i'^' — e^^'^o). Here 
I is the lattice parameter and fsz is the volume of the 
Brillouin zone over which the q integral is performed. 

If the local particle displacements are fully afRne then 
one has U; = DR^, and hence A^ = D(Ri — Rq). Gener- 
ically the displacements will contain a non-affine com- 
ponent and the coarse-grained local deformation tensor 
D can then be defined [H [16] as the one that minimizes 
^JAi — D(Rj — Rq)]^. The minimal value of this quan- 
tity is the non-affinity parameter x- 

To simplify the notation we arrange the Nd rela- 
tive displacement components A^q, where the index 
a = 1, . . . , d labels the spatial components, into an Nd- 
dimensional vector A. We similarly define a vector e 
whose components are the d^ elements of the local de- 
formation tensor, D^^, arranged as a linear array (viz. 
e = (Dii,Di2, . . . ,Did,D2i, ■ ■ ■ ,DNd)), and a matrix R 
of size Nd x d^ with elements Ria,-yY = Sa-fiRi^' — Roy) 
where the Riy and i?o,7' are the components of the lat- 
tice positions R^ and R^, respectively. Below we use the 
notations e and D for the deformation tensor interchange- 



ably, as convenient in the context. 

As explained above, we will define the non-affinity pa- 
rameter x as the residual sum of squares of all the dis- 
placements of the particles in 51 after fitting the best 
affine deformation, measured with respect to the refer- 
ence configuration [8 . The local deformation e is thus 
obtained by minimising the positive definite quantity 
(A — Re)^ with respect to e: 

X = miue (A — Re)^ 
= mine(A*A- A*Re-e*R*A + e*R*Re) (1) 

Here the superscript t denotes the transpose operation. 
The coarse-grained local deformation, i.e. the value of e 
where the minimum is obtained, can then be written as 

e = QA, (2) 

where 

Q = (R*R)-iR*. (3) 
The resulting non-affinity x from ([T]) is 

X = (A- R(R*R)-iR*A)^ = A*PA (4) 

where 

P = I - RQ = I - R(R*R)"iR* (5) 

projects onto the space of A that cannot be expressed 
as an affinc deformation. Note that in arriving at Q we 
have used the fact that P is symmetric, i.e. P* = P and 

p2 = (I - R(R*R)-iR*)2 
= |2 - 2R(R*R)-iR* + 

R(R*R)-i[(R*R)(R*R)-i]R* 
= P. (6) 

As usual this means that all eigenvalues of P are either 
zero or one. 

Having found explicit expressions for e and x '^6 
now proceed to obtain their statistics at low temper- 
ature, where a harmonic approximation will be valid. 
Specifically we consider the canonical distribution of dis- 
placements Uj and momenta at inverse temperature 
/3 = l/fcsT: 



P(p,,u,) = ^cxp[- 13 H{pi,Ui)] 



with the harmonic Hamiltonian 



2to,- 2 



E("» ~ " 



2 

3 J ' 



(7) 



(8) 



where rui is the mass of particle i. The sum in the second 
term in ([s]) runs over all bonds in a harmonic network 
with spring constants K. This is the Hamiltonian for the 
examples we consider in this paper. However, the general 
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expressions that we derive apply directly also to generic 
quadratic Hamiltonians of the form 



p 1 . 

i iaj'j 



(9) 



Here Diaj-y is the dynamical matrix; we have made this 
dimensionless by extracting a factor of l/{f3P). 

Integrating out the momenta from the canonical dis- 
tribution shows that the particle displacements have a 
Gaussian distribution. Their covariances can be ex- 
pressed compactly in terms of the Fourier transform of 
the dynamical matrix \Tl\ 



(10) 



where because of translational invariance the choice of 
reference particle j is arbitrary. This matrix determines 
the variances of the Fourier components according to 



(UqU* > =^"'(q)W(q-q') 



(11) 



where the angled brackets indicate a thermal average. 
The covariances of the displacements are, accordingly, 



J WBZ 



(12) 



For the particle displacements in our coarse-graining vol- 
ume n of interest, — — Uq, we thus find also a 
Gaussian distribution with covariance matrix (A^tjAj^) 
given by 



J VBZ ' 



-iq -Ri 



(13) 



Note that the matrix C defined in this way has the sym- 
metry of the lattice. The thermal average of any observ- 
able ^(A) is then given by. 



(A) 



^ J nrfA,„^(A)exp (-^A*C-^ 



(14) 



with normalization constant — (27r)^''/'^\/det C. In 



the next section we use ( 14 1 to obtain the probability 



distribution functions for x and e. 



III. SINGLE POINT PROBABILITY 
DISTRIBUTIONS 

In this section we derive the single point (local) joint 
probability distribution, P(x,e), for non-affinity x and 



strains e. As before we consider lattices at non-zero but 
low temperatures where a harmonic approximation to 
particle interactions remains valid. To obtain P{x,e), 
we begin with 

which is the characteristic function for the joint probabil- 
ity distribution P{x, e) as measured within Q. Substitut- 
ing the general expressions from Section [nj x = A* PA 
and e = QA, into (fT5| we obtain 



dAia exp 



(16) 



Completing the squares in the argument of the expo- 



nential in ( 16 ) and carrying out the resulting Gaussian 
1 



integrals yields 



<?>(fc, k) = exp |^-^/«*QC(l - 2ifcPC)-^Q*Kj x 

[det(l - 2ifcPCP)]-i/2 (17) 

Setting K ~ and fc = gives the characteristic functions 
of X and e, respectively, as 



<?x(fc) = [det(l - 2ifcPCP)]-i/2 

<Pe(/«) =exp (^-^K^QCQ^H^j 



(18) 
(19) 



Extracting these factors from the joint characteristic 
function shows that it can be written as 

^k,K)^^^{k)^,{K) (20) 

X exp {-ikK*QC{\ - 2ifcPC)-ipCQ*K) 

The term in the second line expresses the fact that x and 
e are generally coupled to each other, rather than vary- 
ing independently. A special case where this does not 
happen occurs when P and C commute. Then one can 
write PCQ* = CPQ*. But this vanishes because from the 
definitions of P and Q one has PQ* = 0. The coupling 
term in ( |2G| then becomes unity and x and e are un- 
correlated. This is the situation we will encounter in the 
one-dimensional example below, when coarse-graining on 
the smallest lengthscale where ft only contains the near- 
est neighbours of particle 0. 

In the case where P and C have a non-zero commutator 
[P,C] = PC — CP, one can put the expansion for small k 



of the coupling factor in ( 20 1 into a form that emphasizes 



the role of this commutator. Specifically, by writing PC — 
CP — [P, C] and exploiting the property PQ* = one finds 



<P{k, k) = •P^{k)<P^{K) exp (-iK*QC [P, C]k + 2i(^ [CP, [P, C]] + [P, C]^) k^ + . . . Q'k 



(21) 
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From the general form (201 of the characteristic func- 



tion, or its expanded version (21 1, we can then obtain the 



desired joint probabihty distribution P(x, e) by inverse 
Fourier transform, either analytically or numerically. 

Before proceeding to apply the above general results to 
two simple example systems, we comment briefly on the 
marginal distributions of x s-nd e whose characteristic 
functions are given in ( 18|19 ) above. From the second of 
these equations, the distribution P(e) of the local strain 
is a zero mean Gaussian distribution with covariance ma- 
trix QCQ*. For the local non-affinity if we call aj the 
eigenvalues of the matrix PCP, then the characteristic 
function ( 18 ) has the explicit form 



1 



(22) 



This shows that x has a generalized chi-square distribu- 
tion P{x)- it is a sum of squares of Gaussian random 
variables, each with zero mean and variance aj. Only 
the nonzero aj contribute here, and there are Nd — cP 
of these. This follows from the fact that P eliminates 
from the space of all relative displacements in ilA the 
d'^-dimensional subspace of affine displacements. 



(a) 
(b) 
(c) 



i = -l Q i = l 



FIG. 1. (a) A portion of a one-dimensional harmonic chain 
showing the neighborhood Q, consisting of particle and its 
two nearest neighbors i = ±1. The affine (b) and the non- 
affine (c) modes are also shown. 



"matrix" in Fourier space is = 2 ^Kl"^ [1 — cos{ql)]] 
in the following we use energy units such that Kf' — 1. 
The displacement covariance matrix ( 13 ) then becomes 



«)(e- 



■iq Xj 



-iqxa 



" 2^ 2^ X '^'^ 1 - cos(qO 

(23) 

which is simply ^ times the identity matrix and so 



A. The one dimensional harmonic chain 

Consider a one-dimensional chain of particles of equal 
mass connected by harmonic springs with spring constant 
K and equilibrium length I as shown in Fig. [T] We choose 
as the coarse-graining neighborhood a central parti- 
cle at Ro = xq = and its two nearest neighbors at 
x±i = ±/. Fluctuating particle positions. Hi = Xi, pro- 
duce displacements Ui — Xi — il and the vector of relative 
displacements is A* — (ui — uq , u_i — uq). The matri- 
ces defined in Section |ll] can be easily evaluated for this 
system and are given by: 



R 



I 



Q = (R'R)-iR* = ( J, -i). 



and 



P = I - R(R*R)-^R* 



1 1 

2 2 



1 1 

2 2 



The two eigenvectors of P corresponding to the eigen- 
values zero and one are (/, — /) and (/, I) respectively. The 
mode with the non-zero eigenvalue corresponds to the 
non-affinc deformation x — (''^i + "-i ^ 2uo)^/2, while 
the one corresponding to the null space of P gives the only 
affine mode e = e = {ui — u-i)/2l of the lattice. The 
affine and the non-affine modes with respect to are 
shown in Fig. [ijb) and (c) respectively. The dynamical 



PCP* = Pp- 



1 1 

2 2 



with eigenvalues a — P/3 ^ and 0, while QCQ* = (e^) — 

rV2. 

The fact that Cij — Pf3^^Sij as found above means 
that the relative particle displacements are uncorrelated. 
This is easy to see intuitively as the potential energy of 
the system is {K/2)^'^^_^{xn+i—Xn—l)'^- Relative dis- 
placements Xn+i — x„ — I of nearest neighbours therefore 
have independent fluctuations, and the relative displace- 
ments Ml— Uq = xi^xq — I and — uq = —{xq — X-i — I) 
in our coarse-graining neighborhood are exactly of this 
form. If we were to enlarge fl, say to include next- 
nearest neighbours, then this would no longer hold as 
e.g. U2 — Uo = X2 — 21 — Xo = {X2 — — /) + {xi — Xq - I) 
is correlated with ui — uq. 



Carrying out the matrix manipulations in (17) after 
specializing to the d = 1 case, we find that the k- 
dependence of the first factor cancels out since [P, C] = 0, 



in line with the general discussion after (21 1. This yields 



the characteristic function for the joint probability dis- 
tribution as the product of the individual characteristic 
functions: 



1 



Vl - 2ika 



(24) 



The joint probability distribution is then obtained by in- 



verse Fourier transforming the characteristic function: 



1 



1 

P(x)P(e) 



exp 



_X_ 
2cr 



(25) 



o 

4 5 



R 



FIG. 2. A typical neighborhood Q around a central particle 
in a triangular lattice, containing the six nearest neighbor 
particles i = 1, . . . , 6. . The reference positions Ro and Ri 
are shown by open circles while the instantaneous positions 
To and Ti are indicated by filled gray circles. 



R and P 



This has a simple form, namely, a product of the chi- 
square distribution of a single Gaussian random variable 
and a Gaussian. We can obtain immediately, for exam- 
ple, the n-th moments of x, (x") = (2cr)"F(n+ i)/F(i), 
which are all finite. In Section |V]we show that one can 
define an external field that couples to x and, for 
a specific value, can cause all the moments to diverge so 
that P(x) crosses over to a distribution with a power-law 
tail. 



B. The two dimensional harmonic triangular net 



The joint probability distribution of local coarse- 
grained strain e and non-affinity x for a two-dimensional 
triangular lattice can be obtained in a similar manner. 
We choose again a nearest neighbor (hexagonal) coarse- 
graining neighborhood fl as shown in Fig. [2] To simplify 
the notation we also assume, without loss of generality, 
Ro = in what follows, and take the lattice constant 
as our length unit so that 1 = 1. Following the lines of 
the 1-d calculation, we begin by obtaining the matrices 
R and P. 



The matrix R is a 12 x 4 matrix encoding the position 
vectors of the 6 neighbors of the particle at the origin 
(see Fig. [2]) . Explicitly one has 



/ Rii 




Rl2 







Rii 



\ 

Rl2 



Rei R&2 

V i?6l i?62/ 



(26) 



Here the index a = 1, 2 indicates the x and y components 
of the lattice positions Ri respectively. 

To find the projection matrix P = I - R(R*R)~iR*, 
we substitute the above form of R into the matrix 
R(R*R)~^R*. One finds that this consists of 6 x 6 blocks, 
each of which is a 2 x 2 diagonal matrix of the form 



1 



{RiiRji + Ri2Rj2) 



I 
1 



R-7 • R7 



1 
1 



The resulting P has 4 zero eigenvalues corresponding 
to the affine transformations. The 8 unit eigenvalues 
correspond to nonaffinc distortions within fl. To iden- 
tify a convenient basis for the non-affine (8-dimensional) 
eigenspace, we choose below the eigenvectors of PCP with 
non-zero eigenvalues. Similarly, a physically meaningful 
basis for the affine (4-dimensional) eigenspace is formed 
by the non-zero eigenvectors of (I — P)C(I — P). 



2. C and PCP 

In order to obtain the statistics of x = A* PA and 
e = QA we need to calculate, as before, the eigenvec- 
tors and eigenvalues of the matrices PCP and QCQ*. As 
discussed above, these a re the matrices determining the 
characteristic functions ( 18|19 1 and hence the marginal 
distributions. We thus require the displacement corre- 
lation matrix C, which in turn is calculated from the 
Fourier-transformed dynamical matrix I'(q). For the 
Hamiltonian Q of a regular harmonic triangular net of 
particles with spring constant K and lattice constant I 
this is 



6 



3-2 cos(gj;) - cos(ig2,) cos(^g.y) 



\/3sin(iga.) sin(^gfy) 
3[1 - cos(igj;) cos(^gj,)] 



(27) 



Here we have again chosen energy units such that KP = 
1. We also use the more intuitive notation = qi and 
qy = q2 for the wavevector components. 

The elements of the real symmetric matrix C are ob- 



tained by evaluating the integral ( 13 ) over the Brillouin 
zone of the triangular lattice. It can be shown, by util- 
ising lattice symmetries, that the integral can be trans- 
formed to one over a rectangular region: 



1 



%/3 



dq^ I dqyV^^{q^,qy) x 
Jo 

e«q-Ro)(e-*q-i^j - e-"i-^o) (28) 



We compute both the real and the imaginary parts 
of this two-dimensional integral numerically, using 256- 
point Gauss-Legendre quadrature. The imaginary parts 
of the elements of C vanish and provide an estimate for 
the accuracy of our numerics. The normalizing volume of 



the unit cell in the reciprocal lattice is vbz 



V3 



. C is a 



6x6 block matrix, with each block of size 2 x 2 as before. 
Of course not all these 36 blocks need to be calculated 
independently, because of the overall symmetry C = C*. 
Using also the additional symmetry relations (see Fig. [2| 

^^1017 — <-^4q47j <-^2q!27 — <-'5a57 ) '-^3a37 — *-^6a67 

c 

Ciq37 ~ ^'4067 J C2a47 = ^^5017: ^^3057 = ^^6027 (29) 

one finds that only 12 blocks of C are distinct. 

With P and C in hand one can construct and diago- 
nalize PCP. This has 12 eigenvalues aj [j = 1, . . . , 12), 
four of which are zero. The 8 non-zero eigenvalues, which 
correspond to the nonaffine distortions within Q shown 
in Fig. [3j are 

/3cti = 2.454 ^ [3(J2 

/3ct3 = 0.482 

/3ct4 = 0.312 = /3cr5 

Pae = 0.283 = I3g7 = /Jcrg (30) 



The structure of the 4-dimensional null space of PCP 
can be understood by looking at the non-zero eigenvec- 
tors of the matrix obtained by the complementary projec- 
tion, viz. (I — P)C(I — P). These eigenvectors correspond 
to the afRne eigendisplacements shown in Fig. |4] 

The independently fluctuating "directions" of the local 
deformation tensor e can be worked out from the covari- 
ance matrix QCQ* of the Gaussian distribution of e (see 
( 19 ). The projections of e onto these directions then have 



familiar forms and map (see below) to the afhne eigendis- 
tortions in Fig.[4j Explicitly we find for these projections: 
(a) volume change (dilation), e„ = ei -I- 64 = Dn -f D22, 
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FIG. 3. Configurations showing the non-zero eigenvectors of 
PCP, which represent non-affine displacements, corresponding 
to eigenvalues in descending order: (a,b) cri — (J2; (c) 0-3; (d,e) 
a4 = 0-5; (f,g,h) (76 = (T7 = erg. 
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FIG. 4. Particle displacements for the non-zero eigenvectors 
of the matrix (I — P)C(I — P), showing afiine eigendisplacements 
(a) dilation (b) uniaxial strain, (c) shear and (d) rotation. 
The reference lattice positions are shown by filled circles. The 
central atom has been deleted for clarity. 



(b) uniaxial strain, e„ = ei — 64 = Dn — D22, (c) shear 
strain, Cs = 62 + 63 = D12 + 1^21, and (d) local rotation, 
a; = 62 — 63 = Di2 — D2i- The associated eigenvalues give 
the relevant compliances for our coarse-graining volume 
O: 



0.261 

0.481 = /3{el) 
0.699. 



(31) 



The statistics of the local deformation tensor P(e) there- 
fore consists of independent Gaussian fluctuations of 
these 4 deformation modes, as illustrated in Fig. [5t 



To finish our discussion of the affine displacements we 
comment briefly on the relation between the eigenvec- 
tors of (I — P)C(I — P), which give the independently 
fluctuating displacement patterns in the affine subspace 
("afRne eigendisplacements"), and the eigenvectors of 
QCQ*, which are the independently fluctuating com- 
ponents of the local deformation tensor ( "eigendistor- 



7 




P(ei 




0.05 0.1 0.15 0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 



FIG. 5. (a) P{x) from the exact calculation (line) com- 
pared with that obtained from molecular dynamics simula- 
tions (points) of a 400 x 400 site harmonic lattice with unit 
particle masses at reduced inverse temperature /3 = 200. The 
system was allowed to equilibrate for 2 x 10^ MD steps with 
a timestep of 10"'^, after which configurations were collected 
for 7 X 10^ MD steps, (b) Plot of P(ei) for the same system 
as in (a). 



tions"). The two matrices are related via 

(I - P)C(I - P) = R(QCQ*)R* (32) 

In the discussion above we treated their eigenvectors on 
the same footing, and indeed each eigendistortion e as an 
eigenvector of QCQ* is related to an affine eigendisplace- 
ment given by A = Re. This works, i.e. Re is indeed an 
eigenvector of ( 32 ) , because R* R is a multiple of the iden- 



tity matrix in the two-dimensional triangular net. This 
simplification will hold in all lattices with sufficiently high 
symmetry. Indeed, one can check that R*R has a. d x d 
block structure where the off-diagonal blocks are zero and 
the diagonal blocks are all equal to ^jR^R*. This di- 
agonal block is a matrix that commutes with the entire 
symmetry group of the lattice so by Schur's lemma 
will be proportional to the identity matrix, unless the 
symmetry group of the lattice is too small. 

Next we look at the statistics of the non-affinity param- 
eter X- The characteristic function ( 18 ) can be written 



out explicitly as in (22 1 where the (t^s are the eigenval- 
ues of PCP. As we saw above, this matrix has eight 
non-zero eigenvalues ui , . . . , erg and four zero eigenvalues 
(79 = . . . = <Ti2 = that do not contribute to (22 1 as they 



correspond to purely affine distortions within VL. The 
eigenvectors associated with the 8 non-affine displace- 
ments are shown in Figj3j Thus, as discussed above in 
general terms, P{x) is the distribution of the sum of the 
squares of Nd — d^ — 8 uncorrelated Gaussian random 
variables, with the variances of these Gaussians given by 
the eigenvalues cti , . . . , erg . A numerical Fourier trans- 
form of ( [22| then gives P{x)- The result is shown in 
Fig. [sja), where we also compare with data from molec- 
ular dynamics simulations; the agreement is evidently 
very good. 

The first and the second moments of x may be obtained 
from successive derivatives of ^x(^) with respect to its ar- 



gument so that (x) 



i-Hd^xik)/dk)^^^ = j:^^,a, 



k=0 



25.426/3-2. 



(33) 



The values for the corresponding quantities obtained 
from our MD simulations of 160, 000 particles and 7, 000 
independent configurations are (6.869 ± 0.005)f3~^ and 
(25.51 ± 0.1)/?-^. These are in excellent agreement with 
our theoretical results; the agreement in all other quan- 
tities shown below is of similar quality. Finally, the n"* 
cumulant of x is given by (l/2)(n — 1)! J^ji'^'^j)" ■ 




0.12 
0.08 
0.04 


-0.04 



0.05 ^ 0.1 



0.15 



FIG. 6. (color-online) (a) Plot of P(x,ei) « P(x)-P(ei) as 
obtained from our calculations neglecting cross correlations 
and (b) the correction term P(x,ei) — P{x)P{e-i), both at 
/3 = 200. Note that the correction term is nonzero though 
small. The results from our simulations for P{Xj^i) s^re in- 
distinguishable from (a). The plots for other components of 
e are similar. 

Having looked at the distributions of local strain e and 
non-affinity x separately, we finally ask about their cor- 
relations. One can verify that, though small, the com- 
mutator [P, C] is non-vanishing, in contrast to the one- 
dimensional harmonic chain case. Indeed, measuring ma- 
trix sizes by the Euclidean norm ||A|| = VTrAA*, we 
obtain for the chosen nearest-neighbour coarse graining 
volume fl 



m\\ = 

/5||[P,C]|| = 
/3||[CP,[P,C]]|| = 



3.875 
0.124 



3.591 X 10" 



(34) 



Neglecting the commutator to first approximation gives 
a joint distribution of x and e that factorizes into P{x) 
and -P(e) without any correlation. The result of this 
calculation for the joint probability P(x, ei) is shown 
in Fig. [6f^a). In actual fact, however, correlations are 
present. This means that the effective compliance of the 
solid depends on the value of x (and vice- versa). The 
effect is quantitatively rather small for our example sys- 
tem, as is clear from Fig. [6f^b) where we plot the first 
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correction to the factorized approximation. More sim- 
ply, the couphng between x ^'^'^ ^ can be assessed by 
looking at cross-correlations like 



(xee*)-(x)(ee*) =2QC[P,C]Q* 



(35) 



This result follows from the fact that the l.h.s. is a third- 
order cumulant because (e) — 0, and so is proporti onal to 
the coefficient of the 0{kK^) term in the expansion ( 21 ) of 
ln<?(fc,K). The commutator [P, C] only appears linearly 
here; higher orders would be needed to express correla- 
tions involving higher moments such as {{x— (x))"ee*)). 



Using the numerically calculated [P, C] in (35) we obtain 
the following third-order cross-correlation of x with the 
eigendistortions: 



ixel) 
ixel) 



{X){el) 
{X){el) 



{xuj^)-{x){co') 
{Xel)~{x){el) 



(36) 
5.108 X 10-^^-'^ 
(37) 



One reads off that non-afhnity is coupled with uniaxial 
strain and shear strain, while local volume strain and 
rotational distortions do not generate non-affinity. (This 
remains true also for correlations involving higher orders 
of X-) The correlations of x with e„ and Cs are positive, so 
that a large affine strain locally is generally accompanied 
by a large non-afRnity x- Conversely, a large value of 
X makes the local region more elastically compliant, i.e. 
typically leads to larger affine strains e. More explicitly, 
the affine displacements (I — P) A conditional on the non- 
affine ones PA can be written as a linear function of 
PA plus Gaussian fluctuations that do not depend on 
PA. Thus one can find the distribution of the affine 
displacements, and hence of e, conditional on x from the 
distribution of the first contribution across all non-affine 
displacements satisfying (PA)^ = x- Because the first 
contribiution is linear in A, one deduces that e is a sum 
of a random contribution proportional to x^^^j ^-nd an 
independent Gaussian contribution. From this it follows, 
for example, that (x"ee*) = (x"+^)Mi + (x")M2 where 
Ml and M2 are two n-independent matrices related to 
C. 

The strength of the coupling between x and e will of 
course depend on the chosen coarse-graining region il, 
and in particular on its radius Rq; we return to this topic 
in Section |VI[ But we believe that the positive correla- 
tion between the strengths of local afRne and non-affine 
deformations is not specific to the harmonic system con- 
sidered here and should hold generally for all solids. 



IV. TWO POINT DISTRIBUTIONS AND 
CORRELATION FUNCTIONS 

We now turn our attention to the spatial correlations 
of X and e. This requires us to consider simultaneously 
the displacement differences in two neighborhoods D, and 
Cl centered on lattice positions Rq and Rq, respectively. 
The vector A is defined as above, with an analogous 



definition for A. The geometry and notation for the 
d = 2 triangular lattice are given in Fig. [Tj the d = 1 
case is straightforward. The local affine strain e — QA 
and non-affinity x — A* PA around Rq are then as before 
for whereas for fl we have the corresponding quantities 
e — QA and x — A* PA. Note that since the reference f2 
is just a translated copy of 57, one has Q = Q and P = P. 

To obtain the joint distribution of x, e, x and e we 
need the joint Gaussian distribution of the displacements 
A and A. These have covariances 



Ci, 



(38) 



While the first two averages are identical, correspond- 
ing to two different but equivalent lattice sites, the third 
quantity encodes the displacement correlations between 
the two different sites. It may be obtained from an ex- 



pression similar to (|13 1 



VBZ 



-iq R(, 



(39) 



Note that Ciaj-y is not symmetric with respect to inter- 
changing Rq and Rq, although the correlation functions 
obtained from it below are, as they must be. 



o o 



o o 



FIG. 7. Typical neighborhoods Q and H around particles 
and in the d = 2 triangular lattice, illustrating the defini- 
tions used for obtaining the two-point correlation functions. 
The labels have the same meanings as in Fig. [2] 

We could now proceed as for the local distribu- 
tion P(x, e) and derive the characteristic function 
(l>{k, K,k, k) of the joint distribution P(x,e, x,e). This 
contains rather too much information to present in a con- 
cise manner, however, so we focus directly on the corre- 
lation functions. The simplest one of these is the strain- 
strain correlator 



(ee*) = (QAA*Q*) = QCQ* 



(40) 



In order to obtain space dependent correlation functions, 
this quantity needs to be evaluated for all choices of Rq 
for fixed Rq (which may again be taken as the origin). 

Alternatively, one may obtain the correlation functions 
in q-space. As a side effect, this avoids the Brillouin zone 
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integration in (39). To see this, write 



r* — / 



/ ^p-,i(q)(eW.-Ro)_i)x 

J VBZ ' 



Now because the reference positions in and are just 
translated copies of each other, one has Rj — Rq = 
Rj — Ro, so that the only Ro-dependence resides in the 
last factor. Defining the Fourier transform of Ciaj-y via 

Cia.jf l^v^2 I t^us rcads 

off 



(e 



17 

-?q-(Rj-RQ) 



!)• 



(42) 



The Fourier transform of the strain correlation functions 
is then simply 



(ee*)(q) = QC(q)Q* 



(43) 



and can be written down in closed form provided the 
dynamical matrix I?(q) for the lattice is known. 

Next we consider the spatial correlation functions of 
the non-affinity, (xx) ~ (x)(x)- K is convenient, at this 
stage, to define the vectors Y — PA and Y — PA with 
components yj and ijj where j = 1, . . . , Nd. Thus x = 

X^jlfi y| ^^"i XX = Tl!ij=iyiyj ■ Hence the correlation 
between x and x is given by 



Nd 



Nd 



Nd 



= 2Y,{y,y,r = 2E(PCP)^ 
= 2Tr(PCP)(PCP)* (44) 



using Wick's theorem. If in line with our earlier notation 
we use a'j to denote the Nd eigenvalues of the matrix 
(PCP)(PCP)*, then the last expression can be simplified 
to 



Nd 



(Xx)-(x)' = 25]a 



(45) 



Note that if PCP itself happens to be symmetric, then the 
(t| can be obtained as the squares of the eigenvalues of 
this matrix. As for the real-space strain correlator, one 
has to evaluate ( 45 1 for differt choices of Rq to obtain 



spatial profiles of the non-affinity correlator. 

Finally one could ask about spatial cross-correlations 
like (xee*) — (x)(ee*). We do not pursue this here: as 
we saw above, these correlations are already rather weak 
(at least for coarse graining across nearest neighbours) 
locally, i.e. when Rq = Rq- 



A. The one dimensional harmonic chain 

We now apply the above framework to the one- 
dimensional harmonic chain introduced in Sec. IIII Al We 
choose as the first reference location xq = and as the 
second xq = ml. Coarse graining will be across the near- 
est neighbours x±i in Q and x±i = (to ± 1)^ in Q. The 
matrix Cjk is then given by. 



2-/' dq_ F(g) 
27r// V{q) 



where 



F(g) = (e^«^^ - l)(e-*« = 



— e 



iq Xq \ 



(46) 



(47) 



-iqkl 



1 



There are now two possibilities, li j — k then, bearing in 
mind that 'D{q) = 2§[l - cos{ql)], one has F{q)/'D{q) = 
^--ig-^igso so that Cjk = except when xq = 0. Oth- 
erwise, if J = -k, then F{q)/D{q) = -;3-ie'9(j'-*o)_so 
that now Cjk = unless = jl- For all other cases, fife 
vanishes identically. Summarizing, C equals C from (23) 
when X{) = = xq] for X{) — I it is given by 



C = 



-1 




(48) 



while for xq = —I one obtains the transpose. For all 
larger distances |a;o| > 2Z, C = 0, indicating that x 
and e are uncorrelated beyond nearest neighbours. The 
intuition here is as discussed after (23), namely that 



the relative particle displacements of all nearest neigh- 
bour pairs fiuctuate independently from each other. Ac- 



cordingly, the single nonzero entry in ( 48 ) conies from 
the correlation of the displacements ui = xi — xq and 
M_i = X-i — Xq — Xq — xi, and SO is simply the negative 
variance of mi. 

To obtain the correlation functions we need the ma- 
trices PCP* and QCQ*, where we can focus directly on 
the only non-zero correlations at xq = ±1. Even though 
C is then not symmetric, PCP* is, and has eigenval- 
ues a = ^P/3~^ and 0. On the other hand, the scalar 
QCQ* equals Using (451, the correlation functions, 

(x(0)x('™0) ~ (x)^ and {e(0)e{ml)), can then be written 
as follows: 



(X(0)x(™0)-(X)' = 2^2 = 2/4/3-2, 
= 2a2 = i/4r^ 



= 0, 



(6(0)6(toO) = 



TO = 

m = ±1 
|to| > 2 

TO = 

m = ±1 
ItoI > 2 
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As expected the correlation functions, being symmetric, 
depend only on the magnitude and not on the direction 
of Xq = ml. 

Summarizing, the correlation functions for x and e in 
d = 1 are always short ranged, vanishing identically be- 
yond the nearest neighbor. Correlation functions in the 
two dimensional lattice have much more structure as we 
will see next. 



B. The two dimensional harmonic triangular net 

The calculation of the correlation function for the two 
dimensional triangular lattice follows along lines similar 
to that of the single point distribution functions, except 
that results now have to be obtained for each lattice po- 
sition Rq. For the x correlations, we calculate the trace 
of (PCP)(PCP)* as the sum of the 12 eigenvalues ctJ and 

(x)^, where R = Rn — Rn 



use (45) to find (x(O)x(R)) — yXTi wnere rv, = no — no. 
The results obtained using this calculation are compared 
with simulations in Fig. [Sj showing that this function is 
isotropic and decays within about 2 — 3 lattice spacings, 
similar to the results obtained in Refs. ^Tl I22| . 




FIG. 8. (a) Surface plot of /3^((x(0)x(R.)) - (x)^) over the 
two dimensional triangular lattice from our exact calcula- 
tion, (b) Quantitative comparison for the orientation aver- 
aged l3^{{x{0)x{R)) - ixf), where ii = |Ro - Ro|. Open 
circles show the values obtained from our exact calculation 
and filled circles those from simulations of the 100 x 100 par- 
ticle system. The solid lines are guides to the eye. 

Similarly the spatial correlation function for the strain 
may be obtained by evaluating (ee*) = (QAA*Q*) = 
QCQ* for a range of spatial separations R. The results 
are shown in Fig.|9] The correlation functions of and uj 
decay rapidly to zero and are nearly isotropic. We take 
advantage of this approximate symmetry by averaging 
over all pairs (Ro,Ro) related by symmetry to produce 
angle-averaged correlation functions that are functions of 
R = IRq — Ro| alone. Results for these functions are com- 
pared with those obtained from simulations in Fig. |10[ a). 
Correlations involving the uniaxial (e„) and shear (e^) 
strains exhibit a pronounced 4-fold anisotropy at large 
distances, with prominent lobes at 0, tt, ±7r/2, ±7r/4 and 
±37r/4. Furthermore, (e„(0)e„(R)) and (es(O)es(R)) ap- 
pear to be rotated by tt/A with respect to one another (see 
Fig. [9]jb,c)) for large |R|. At small distances, of course. 




FIG. 9. (color-online) Density plot of the correlation func- 
tions as obtained from our calculations showing the shape 
of the correlation functions in the two dimensional plane 
(a) /3(e„(0)e„(R)), (b) /3(e„(0)e„(R)), (c) /?{e.(0)e.(R)), (d) 
/3{cj(0)a;(R)). Note that the volume and rotational correla- 
tions are nearly isotropic while the uniaxial and shear strain 
correlations show four-fold anisotropy. The colors vary from 
dark blue (—0.1) to white (0.2) for all the plots. To keep a 
uniform scale for all graphs we have cut off large values at the 
origin where necessary. 



this correspondence is not exactly satisfied because the 
underlying triangular lattice does not have this 7r/4 ro- 
tational symmetry. 

An identical 7r/4 anisotropy is also observed in momen- 
tum space (Fig. 11 ) where one can obtain closed form ex- 
pressions for the correlation functions, in particular in the 
q — > limit as we show below. We begin by writing down 
explicit expressions for the strain correlation functions in 
component form. With the abbreviation E = QC(q)Q* 
one has 



(e?)(q) = 


Eiiii - 


- E2222 


+ 2E1122 


(e^>(q) = 


Eiiii - 


- E2222 


— 2E1122 


(e?)(q) = 


E1212 - 


- E2121 


+ 2E1221 


^'>(q) = 


E1212 - 


" E2121 


— 2E1221 



(49) 



The notation used here is the same as in (43 1, so that e.g. 
(e^)(q) is the Fourier transform of the strain correlator 



(ei,(Ro)e„(Ro)) 



(q))e 



,iq-(Ro-R.o) 

After substituting in for C(q) from ( |42[ ) one gets, using 
also that for the triangular lattice R*R is three times the 
identity matrix. 



1 ^ 



(50) 
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FIG. 10. (a) Quantitative comparison for the (orientation 
averaged) volume and rotational deformation correlations 
/3(e„(0)e„(_R)) (squares) and /3{aj(0)a;(i?)) (circles). The open 
symbols are from our exact calculations and the filled ones 
are from our simulations of the 100 x 100 particle system. 
The solid lines are guides to the eye. The agreement in the 
numerical values for the correlations for other components of 
the distortion tensor are similar, (b) The large length scale 
behaviour of /3{e„(0, 0)eu(_R/^/2, ^/v^)) (open circles) and 
l3{es{0,0)es{0, R)} (filled circles) showing the slow decay along 
these directions. The dotted line shows the ~ behaviour 
for comparison. 



with 

fi = cos[q- (Ri 
f/ = sin[q • (R„ 



- R-ri) 



cos(q • R;) — cos(q • Rj) 
~ sin(q • Ri) + sin(q • Rj 



(51) 
+ 1 



Noting that the imaginary c ont ributions to E sum to zero 



and expanding /]j in Eqn. ( 50 ) yields for small q 



(55) as 



(e?(q)> 



64 



9(9^ 



qy\^fqx- qy ^ ^ 



V2 



V2 



(57) 



This evidently maps to ([54]) under a rotation by 7r/4, 
as claimed. We note that the large distance (small q) 
anisotropics of (e^(q)) and (e^(q)) are also consistent 
with a mean field, continuum theory, calculation shown 
in detail in Ref. Interestingly, that calculation re- 
stricted attention to fluctuating strain fields that satisfy 
force balance, so one concludes that it is indeed these con- 
figurations that dominate the scaling for large distances. 
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2 

qy 
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Qx 



1 

g„'<Zyp-l(q) (52) 



id. 



To 



where we have used again that RiaR 
examine the leading order behaviour at small q for the 
strain correlators, we expand I?^^(q) about q = in (52) 
and substitute into ( 49 1 to obtain 



(e^Xq) 

(eS)(q) 

(e^)(q) 
(^')(q) 



64 



9 ? 

qiqy 



9 iql 
64 

^ (ql + ql 



9 9 

qUv 



(53) 
(54) 

(55) 
(56) 



consistent with the results plotted in Fig. [TT] The fact 
that the correlators of uniaxial strain e„ and shear strain 
Cs are related by a rotation becomes clearer if one rewrites 



FIG. 11. (color-online) Plot of the correlation functions in 
q- space obtained from (421. (a) (e^)(q), (b) {e^)(q), (c) 



(e3)(q), (d) {uj )(q). The key to the colour values are given 
beside each plot. Note the 4-fold symmetry of the uniaxial 
and shear correlations similar to the real space plots shown 
in Fig.[9l 

The q space structure of the uniaxial and the shear 
correlation functions as implied by (561, Figs. [TT|b) and 
(c), and further supported by continuum theory |22j . 
shows that these correlation functions have singularities 
at g = 0, with the second terms in (54) and (57) van- 



ishing along specific directions (qx = 0,qy = in e„ 
and = iqy in e^). These singularities lead to slow 
(~ 1/R^) decay of the correlation functions in real space 
(see Fig. [To]^b)), with the prefactor alternating in sign 
according to cos(4(?) or sin(40) with the polar angle of 
R. 

Before we end this section, we remark that the strain 
correlation functions, by linear response, are proportional 
to the strain field produced by a point, delta function 
stress at the origin. Since the large R (or small q) struc- 
ture of the correlation functions is insensitive to crystal 
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symmetry, it is no surprise that similar quadrupolar dis- 
placement patterns |30j have been observed in associa- 
tion with local rearrangements in amorphous materials 
known as shear transformation zones (STZ) "f, '5', in 
both experiments ,3 Ij and computer simulations .32i 
of granular or glassy materials under shear. 



LINEAR RESPONSE AND THE 
NON-AFFINE FIELD 



The form of the characteristic function ( 15 ) offers a 



simple way to calculate the response of the system to 
uniform fields conjugate to x ^'^^ ^- We analytically 
continue k and k to the complex plane by replacing 
K — > K—i S and k — )■ k—i h^. Here the vector S, once re- 
arranged into a symmetric tensor, is the stress [5] , made 
dimensionless by multiplying by the inverse temperature 
/3 and the size of a suitable local volume of the order of 
Rf^. On the other hand is a new field, conjugate to 
X- The introduction of a (small) stress merely shifts (e) 
away from zero to a value proportional to S (Hooke's 
law), i.e. (ei)s = (eiej)x;=oSj. The proportionality con- 
stant here is the zero field compliance calculated earlier. 
Furthermore, because of the coupling between e and x, 
external stress (mainly shear and uniaxial, see Section III 
and Eq. (36 1) will change P{x) for lattices where [P, CJ 
and higher commutators are non-zero. A straightforward 
calculation, introducing S in ( |21[ ) and Taylor expanding, 
shows that for small values of S, 

(x)i: = (x)5:=o + S*QC[P,C]Q*S, 

which always increases (x) for the d — 2 triangular lat- 
tice; higher moments of x are similarly affected. At large 
S, of course, perturbations would become so large that 
the effects of anharmonicities that we have not modelled 
would become apparent. 

The effect of is also intriguing. Below we illustrate 
this explicitly for the one-dimensional case, though simi- 
lar results should hold in any dimension and for particles 
with arbitrary interactions. The joint characteristic func- 
tion for X a-nd e for the d = 1 chain with included is 



^(fc, k) 



1 - 2 cr /;,v 



■ cxp 



1-2(7 ifc-2cr/i^ V 2 



1 



2\ 2 



(58) 



Note that we have multiplied <?(fc, k) by the fac- 
tor y^l — 2a hy, to ensure normalization, i.e. 
limfe_i.o.K-i.o^(fcj k) = 1- Now we can obtain P(e), 
P{x) and P{x, e) by inspection as 



P{x,e)^Pix)P{e) 



l-2ah 



2na 



-X ^^^exp 



{l-2ah^)x 



2a 



exp 



2(e2 



(59) 



As — ^ l/(2cr), P{x) becomes proportional to x 
and all the moments 



-1/2 



(x"> 



2a 



l-2aK 



^(r^+^) 

r(i) 



diverge as [/i^ — l/(2a)]~". Spatial correlations of non- 
affinity X also become long ranged in this limit and the 
system becomes disordered. Displacements acquire a 
non-affine character over all length scales, as evidenced 
from the decrease of the amplitude of the structure fac- 
tor for a finite one-dimensional chain of 100 particles, see 
Fig.pta) and (b). 




(b) 



FIG. 12. (a) The scaled, two point X'X correlation function 
for = 0, 0.4 and 0.49 from Monte Carlo simulations of a 50 
particle harmonic chain at /3 = 100, showing a large increase 
in the correlation length as the critical field is approached, 
(b) The first peak in the structure factor S{q) — (pqP-q) for 
= (solid line), 0.4 (dashed line) and 0.49 (dash-dotted 
line) for the same system as in (a). Note that the amount of 
structure is reduced at constant temperature. 

All of this suggests the presence of a critical line in the 
(/i^,/3) plane beyond which the system becomes globally 
non-affine so that x defined by coarse graining over local 
neighborhoods f2 is always infinite - a "maximally non- 
affine" solid (the presence of anharmonicities would, in 
practice, be expected to limit x to a finite value). For 
the d — 1 lattice, this transition is identical to the cel- 
ebrated Peierls transition as can be deduced from 
the nature of the non-affine mode shown in Fig. [ijc) . In 
higher dimensions, the transition appears for values of 
equal to a critical ft.* close to half the reciprocal of 
the largest eigenvalue of PCP*. Indeed, one finds easily 
by writing down the general analog of ( 58 ) , taking a log 
and differentiating w.r.t. ik that 



(X) - Tr(l 



(60) 



which diverges at /i* = l/{2niaxj aj) if the aj arc the 
eigenvalues of PCP* as before. In the triangular net, the 
largest eigenvalues are CTi = CT2 and the displacement pat- 
terns whose amplitudes would grow at the transition are 
shown in Fig. [3](a) and (b). Locally, they correspond to 
an almost uniform translation of the neighbouring parti- 
cles relative to the central particle, though it is difficult 
to visualise what global configuration is finally produced. 
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We imagine that this leads to destruction of the lattice 
structure and eventual amorphisation. Of course, as the 
transition is approached, we expect the identity of the 
neighbourhoods, SI, itself to become ill-defined due to 
this loss of crystalline order, making many of our results 
invalid in that limit. 

In the harmonic lattice, the transition discussed above 
is hidden because the only physically realisable value of 



the non-affine field 
finite temperature /3 



- lies on the critical line at in- 
0. Nevertheless, there may be 
systems where the critical line cuts the = axis at a 
non-zero value of /3 (i.e. at finite temperature). In such 
a system one may obtain a physical transition from an 
affine to a maximally non-affine state as the temperature 
is increased or a stiffness parameter is reduced producing 
a going over to a "glass spinodal" [34] . We speculate that 
this may also happen at or near the yield point of a solid 
under external load where the bonds between atoms be- 
come weak due to strong anharmonicities. However, for 
such cases, the simple linear response calculations pre- 
sented above becomes invalid and the distributions of x 
and e become strongly coupled through the strain and 
non- affinity dependence of the dynamical matrix 2?(q). 
Implications of this transition for the mechanical and 
phase behaviour of solids in 1, 2 and 3 dimensions are 
being worked out and will be published elsewhere. 



Firstly, x depends on the number of sites within fl, so 
we need to consider the intensive quantity (x) /Rq where 
i?o is the radius of the reference volume fl. Furthermore, 
since x depends on the fluctuations of the displacements 
u, then in d = 1 this quantity itself should diverge as 
~ Rn and in d = 2 as log(i?n) [1^. In Fig. [l3| we 
show plots for this normalized non-affinity in d = 1 and 
d = 2, which confirm these expectations. In higher di- 



ix) 




FIG. 13. Plot of {x)/Rn vs Rn for d = 1 (filled circles) and 
d — 2 (open circles). Note that while the intensive (per unit 
coarse-graining volume) x increases linearly with Rq, the size 
of the reference volume SI, in d = 1, it has a much slower 
logarithmic increase in d = 2. 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have shown that coarse graining of the 
microscopic displacements of crystalline solid at non-zero 
temperatures generates non-affine as well as locally affine 
distortions. The procedure effectively amounts to inte- 
grating out phonon modes with wavelengths comparable 
to or smaller than the coarse-graining length. We have 
obtained the probability distributions for the non-affine 
parameter and the affine distortions of a harmonic lattice 
at non-zero temperatures. We have also obtained the spa- 
tial correlations of the local distortions and non-affinity. 
While (x(O)x(R)) decays exponentially, the correlation 
functions corresponding to the different elements of the 
distortion tensor decay differently. Volume and rotation 
correlations, on the one hand, are short ranged; uniaxial 
and shear strain components, on the other hand, decay 
as a 1/7?^ power law with prefactor depending on the 
polar angle as sin(46') or cos (46*), respectively. The angu- 
lar dependences of the slow decay for uniaxial strain and 
shear are rotated by 7r/4 with respect to each other. We 
noted that this is consistent with an earlier continuum 
calculation imposing the mechanical stability condition 
that the stress is divergence free in the fluctuating strain 
field. Finally, we have shown that it is possible to induce 
a transition from a solid with a finite average non-affinity 
(x) to one where all moments (x") diverge, by tuning a 
non-affine field h^. 

How do our results vary with the size of the reference 
volume n used to define the coarse-grained quantities? 



mensions the intensive variable (x) /Rfi should for large 
i?o approach a constant (proportional to temperature, as 
for d = 1 and 2). Note that the probability distribution 
of strains, P(x,ej), is similarly fl dependent and so our 
calculation automatically incorporates exact finite size 
scaling of the elastic compliances. Approximate finite 
size scaling results based on continuum elasticity theory 
have been used to obtain elastic constants of colloidal 
solids from video microscopy data [201 121] • Our results 
may offer a better way to analyze such data. 

Secondly, as Rn increases, x and the distortion D may 
get more and more coupled. This is best illustrated, for 
the d — I harmonic chain, by computing the norms of 
the commutators [P, C] and [CP, [P, C]] which we obtain 
as shown in Table |lj 

For very large Rq more and more terms are needed to 
get good convergence of the Taylor expansion for P(x, ^j) 
in powers of k, viz. Eq. (21), for given k. Accordingly 
X becomes inextricably linked with the strain, a phe- 
nomenon related to the fluctuation-driven instability of 
ordered solids in one dimension [lOj. This effect should 
be weaker in higher dimensions, though a full study of the 
influence of dimensionality on non- affinity for a variety 
of lattices and interactions needs further work. 

Our calculations may be easily extended to other lat- 
tices and to higher dimensions without much difficulty, 
requiring at most a calculation of the relevant dynam- 
ical matrix. Similarly, local x E^nd displacement distri- 
butions can be obtained for crystals with isolated, point 
or line defects once the appropriate Hessians of the local 
potential energy at defects are evaluated. The effect of 
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l|C|| 


|[P,C]|| 


ll[CP,[P,C]]|| 


1 


1.414 


0.000 


0.000 


2 


3.741 


0.282 


0.126 


3 


7.211 


0.654 


0.534 


4 


11.832 


1.148 


1.516 


5 


17.606 


1.766 


3.448 


6 


24.535 


2.507 


6.803 


7 


32.619 


3.371 


12.145 


8 


41.856 


4.358 


20.132 


9 


52.249 


5.469 


31.520 


10 


63.796 


6.704 


47.154 



TABLE I. Numerical values of the norm of the successive 
commutators of P and C which contribute to the cross cor- 
relation of non-afBnity and strain (see section Note that 
the commutators increase with the size Rn of the reference 
volume. 

external stress on x is another interesting problem which 
may be addressed immediately in the limit of negligible 
anharmonicity; the effect of anharmonic terms could be 
included perturbatively. Finally the effect of disorder can 
be incorporated [12j at non-zero temperatures. 

Our results may also be used to construct new simula- 
tion strategies for investigating the mechanical behaviour 
of solids under external stress. For instance, particle 
moves may be designed which change x without influ- 
encing the local distortion within by projecting particle 
displacements along eigendirections of P. 

Such calculations will be particularly useful for 
glasses jT,, where local potential energy Hessians may 
be used to define an equivalent harmonic lattice at ev- 



ery instance of time with x being calculated dynami- 
cally from the reference configurations at the previous 
timestep. Calculations similar to ours will be useful to 
understand the properties of shear transformation zones 
(STZ) [29j , defined ^|8| as regions with a large value of x; 
these are the dominant entities responsible for mechani- 
cal deformation of glasses. STZ are thermally generated 
and respond to external stress by rearrangements of lo- 
cal particle positions. Similar localised non-affine excita- 
tions have also been observed in anharmonic, crystalline 
solids 25 where they have been identified as droplet fluc- 
tuations from nearby glassy and liquid-like minima of the 
free energy. Constrained simulations like those outlined 
above may help in identifying the role of x in the pro- 
cesses involved in complex phenomena like anelasticity, 
yielding and melting. The role of the non-affine field hy- in 
influencing glass transition and the mechanical behaviour 
of solids both crystalline and amorphous is another di- 
rection that we intend to investigate in the future. 
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